* Replication of Figure 8, panel b, using Golder marginal effects plot
* You will need to specify the directory in which the repliaction data are stored in the second line of code 

clear
cd "[Directory with replication data here]"

version 11.0
#delimit ;
log using "sim8b.log", replace;
set more off;

use "fordham_wp2025.dta", clear ;

logit multilat time black interactb ag interacta lynchings interactl cotton textapp, cluster(icpsr);

preserve;

set seed 339487731;

drawnorm SN_b1-SN_b10, n(10000) means(e(b)) cov(e(V)) clear;

postutil clear;
postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi 
            using "sim8b.dta", replace;
            
noisily display "start";

local a=0 ;
while `a' <= 70 { ;

{;

scalar h_time=0;
scalar h_cotton=0.26;
scalar h_textapp=4.7;
scalar h_black=20;
scalar h_ag=20;
scalar h_lynchings=30;
scalar h_constant=1;

    generate x_betahat0  = SN_b1*h_time 
                            + SN_b2*h_black   
                            + SN_b3*h_time*h_black 
                            + SN_b4*`a' 
                            + SN_b5*`a'*h_time
                            + SN_b6*h_lynchings
                            + SN_b7*h_lynchings*h_time
                            + SN_b8*h_cotton 
                            + SN_b9*h_textapp
                            + SN_b10*h_constant;
                            
    generate x_betahat1  = SN_b1*h_time 
                            + SN_b2*h_black  
                            + SN_b3*h_black*(h_time+17)
                            + SN_b4*`a' 
                            + SN_b5*`a'*(h_time+17) 
                            + SN_b6*h_lynchings
                            + SN_b7*h_lynchings*(h_time+17)
                            + SN_b8*h_cotton 
                            + SN_b9*h_textapp
                            + SN_b10*h_constant;
    
    gen prob0 =normal(x_betahat0);
    gen prob1=normal(x_betahat1);
    gen diff=prob1-prob0;
    
    egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
    
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 diffhat;
    
    local a=`a'+ 1;
    
    display "." _c;
    
} ;

display "";

postclose mypost;

restore;

merge using "sim8b.dta";

gen yline=0;

gen MV = _n-1;

replace  MV=. if _n>70;

graph twoway hist ag, percent color(gs14) yaxis(2)
        ||  line diff_hat MV, clwidth(medium) clcolor(blue) clcolor(black)
        ||  line diff_lo MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line diff_hi MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line yline MV,  clwidth(thin) clcolor(black) clpattern(solid)
        ||  ,   
            xlabel(0 20 40 60 70, nogrid labsize(3)) 
            ylabel(-1 -0.8 -0.6 -0.4 -0.2 0 0.2, axis(1) nogrid labsize(3))
            ylabel(0 5 10 15 20 25 30 35 40, axis(2) nogrid labsize(3))
            yscale(noline alt) yscale(noline alt axis(2)) xscale(noline) legend(off)
            yline(0, lcolor(black))
            yline(-1 -0.8 -0.6 -0.4 -0.2 0 0.2, lcolor(white))
            xtitle("Percentage of District Workforce in Agriculture", size(3))
            ytitle("Marginal Effect of Moving from 1945 to 1962", size(3))
            ytitle("Percent of Observations" , axis(2) size(3))
            xsca(titlegap(4)) ysca(titlegap(4)) ysca(axis(2) titlegap(4))
            scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));
            
log close;
